Velocity Fluctuations in Dynamical Fracture: the Role of Microcracks 
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We address the velocity fluctuations of fastly moving cracks in stressed materials. One possible 
mechanism for such fluctuations is the interaction of the main crack with micro cracks (irrespective 
whether these are existing material defects or they form during the crack evolution). We analyze 
carefully the dynamics (in 2 space dimensions) of one macro and one micro crack, and demonstrate 
that their interaction results in a large and rapid velocity fluctuation, in qualitative correspondence 
with typical velocity fluctuations observed in experiments. In developing the theory of the dynamical 
interaction we invoke an approximation that affords a reduction in mathematical complexity to a 
simple set of ordinary differential equations for the positions of the cracks tips; we propose that this 
kind of approximation has a range of usefulness that exceeds the present context. 
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I. INTRODUCTION 

Classical linear elasticity fracture mechanics provides 
clear cut predictions for the dynamical evolution of cracks 
in stressed materials. Under pure mode I loading a crack 
is expected to remain straight, and to exhibit a tip ve- 
locity that increases monotonically towards the Rayleigh 
wave speed cr. Reality shows that this is but a pipe 
dream. When the crack velocity exceeds a finite fraction 
of cr the velocity of typical cracks exhibits wild fluc- 
tuations, the crack surfaces lose their smoothness and 
the mean velocity never asymptotes towards cr. The 
fundamental understanding of the discrepancy between 
the prediction of the classical theory and experiments 
remains an open problem of considerable interest and 
importance. 

A number of studies 

0,11111111111111 point towards 
a close correspondence between the onset of velocity fluc- 
tuations and the appearance of secondary damage like 
micro cracks (appearing ahead of the tip), microscopic 
side branches etc. Conical markings which are observed 
on crack surfaces offer a good indication that micro cracks 
exist before the arrival of the crack, although it is not 
determined whether the former stem from material im- 
perfections or from stress instabilities. The fact that the 
density of conical markings increases during the crack 
evolution [f| suggests that the level of stress is responsi- 
ble in some way for the activation of the micro cracks. 
The aim of this paper is to explore the connection be- 
tween velocity fluctuations and the putative existence of 
micro cracks ahead of the crack tip. To this aim will 
study the dynamical interaction between a macro crack 
and a micro crack and focus on the velocity of the tip of 
former under the influence of the latter. 

To actually solve exactly the dynamical equations for 
the displacement field with boundary condition on both 
macro and micro cracks up to coalescence is a very taxing 
quest. Building upon experience in the field we will pro- 
pose here an approximate methodology that will allow us 
writing down ordinary differential equations for the posi- 
tions of the tips of both macro and micro cracks. While 



sensible, the approximate methodology is not established 
in a controlled fashion, requiring therefore simulational 
support. Indeed, we will offer in this paper lattice simu- 
lations to back the analytic considerations. We will show 
that the correspondence is excellent. 

In Sect. [n]wc introduce the problem at hand, being an 
infinite 2-dimensional stressed material with one macro 
crack and one colinear micro crack of length I. In Sect. 
IIIII we describe the approximate method of solution, mo- 
tivating it by the exactly soluble cases of straight and 
bifurcating cracks. The Section culminates with approx- 
imate equations of motion for the tips of the macro and 
micro cracks. In Sect. II VI we describe the solution of the 
model problem, stressing the velocity of the tip of the 
macro crack. We show that the net result of the inter- 
action is a rapid and large up and down fluctuation in 
this velocity, in correspondence with the observed fluctu- 
ations in dynamical crack propagation. Sect. Ivl provides 
a simulational support to the approximate theory; by 
performing lattice simulations we study the same model 
problem and compare the results. The close correspon- 
dence between approximate theory and simulations lends 
support to the former. Sect. IVII offers a summary and 
conclusions. 



II. THE PROBLEM 

The problem that we want to consider is sketched in 
Fig. n We consider a macro crack and a micro crack 
that at a given time extend along the intervals [— L, 0] 
and [^_,^ + ] respectively. The distance between them is 
given by A. The length of the micro crack is £ = £ + — £-. 
We expect on physical grounds that the micro cracks in 
typical materials are at most of the size of the process 
zone, and therefore we always consider the limit £/L — > 0. 
The aim of the calculation is to determine the simulta- 
neous motion of the three crack tips (the macro crack tip, 
the inner and outer tips of the micro crack) as a function 
of time. In full generality this entails the general solution 
of the field equations for an arbitrary motion, specifically 
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FIG. 1: The geometric configuration of the model problem. 
The macro crack and the micro crack extend along the inter- 
vals [— L, 0] and respectively with Lji 3> 1. 



the determination of the dynamic stress intensity factors 
at the cracks tips and then to apply a fracture criterion 
to obtain the actual dynamics. We cannot offer an ex- 
act solution to this problem. Instead, we will introduce 
an approximate method that provides analytic insight to 
the problem. 



III. APPROXIMATE METHOD OF SOLUTION 

To motivate our approximate methodology we will re- 
call some exact classical results obtained for ideal (mode 
I) straight dynamical cracks and more recent results per- 
taining to (mode III) bifurcating cracks. 

A. Motivation I: Ideal Straight Cracks 

As said above, linear elasticity fracture mechanics pro- 
vides exact solutions for straight cracks under mode I 
loading. The stress field <7ij(r, i), measured in polar 
coordinates relative to the tip, is expected to have a uni- 
versal form in the vicinity of the crack tip 9] , 



a ij {r,6,t)=K 1 {v{t),L(t)) 



^■(0, «(t)) 



(1) 



Here v(t) and L(t) are the time-dependent crack veloc- 
ity and length respectively, (9, v) are known universal 
functions Q, and Kj(v, L) is the "stress intensity fac- 
tor" which is predicted to depend on the instantaneous 
crack velocity and length only. For notational simplicity 
we drop the t dependence. At each moment in time the 
velocity is expected (for plane stress conditions) to be 
determined by the energy balance equation 



r(v) = -A 1 (v)K?(v ) L) 



(2) 



The LHS here is the fracture energy and the RHS is the 
energy release rate into the crack tip region, resulting 
from a path integral over the total energy flux. E is 
Young's modulus and A Y (v) is a mode I universal func- 
tion. The exact result that we refer to is the decompo- 
sition of the dynamic stress intensity factor K I (v, L) for 
a semi-infinite crack under time independent loading, in 
the form [gl fiol 



where kj(v) is a universal function of v and K^(L) is 
the stress intensity factor of a static crack of length L 
under the same loading (when L is large enough to be 
considered as semi- infinite) . This important result is the 
basis of the classical theory of straight crack motion. The 
calculation of the static stress intensity factor is a much 
easier task than the evaluation of its dynamical counter- 
part, since it requires solutions of bi-Laplace equations 
with boundary conditions. Rewriting Eq. J2J in the light 
of this result one obtains, 



T(v) = ±A 1 (v)[k 1 (v)K°(L)f 



(4) 



A further serendipitous simplification arises from numer- 
ical evaluations of the combination A 1 (v)k'^(v), showing 
that it is well approximated by 



A T (v)k*(v) re l-v/c R . 



(5) 



This approximation leads to an ordinary differential 
equation for the crack length. If one asserts that the 
fracture energy T is u-independent this differential equa- 
tion becomes explicit, 



dL(t) 
dt 



cr 



ET 



(6) 



It should be stressed again that the decomposition prop- 
erty of the dynamic stress intensity factor (Eq. |J3J) is 
essential in deriving this basic equation. This equation 
predicts a monotonic increase in the tip velocity asymp- 
toting towards cr. There is no crack motion as long as 
the stress intensity factor does not exceed a material de- 
pendent threshold 



K*(L{t)) < 



no crack motion . 



(7) 



For the sake of illustration we show in Fig. the so- 
lution of this equation for the simple case K^[L(t)) = 
a °° y/irL(t)/2, where a°° is the load at infinity. 

B. Motivation II: Bifurcating Cracks 

The bifurcation of fast cracks is observed in many ex- 
periments, and understanding it theoretically is a prob- 
lem of some importance in the theory of fracture. One 
interesting problem that was addressed recently [T^ in 
this context is the determination of the stress intensity 
factors at the tips of symmetrically branched cracks in 
terms of the stress intensity factor prior to branching. 
Ref. 01 presented a solution to this problem for mode 
III (anti-plane) conditions. In this solution the stress in- 
tensity factor K' at the tips of two symmetric branches 
emerging from a macro crack at a velocity v' and creating 
an angle Xn relative to the macro crack line, was given 
in the form 



K l {v,L)=k l {v)K°{L) 



(3) 



K' = y/1 - v'/c a H 33 (X,v'/c s ) K 



(8) 
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FIG. 2: The predicted velocity increase as a function of nor- 
malized time (where Lq is the crack length at initiation) for 
a mode I crack under uniform constant load at infinity, Eq. 



Here y/l — v'/c s is the mode III universal function, 
whose mode I counterpart is ^(v), c s is the shear wave 
speed, Kq is the stress intensity factor of the macro crack 
prior to branching and H 33 {\, v' /c s ) carries the informa- 
tion regarding the dynamic interaction. Note that the 
macro crack velocity v is absent here as in the branching 
scenario adopted in ^1 the macro crack stops suddenly 
before branching and a static stress distribution is estab- 
lished behind a wave front travelling at the characteristic 
wave speed c s . If the decomposition in Eq. (|3J is of some 
generality then we expect the main interaction effect to 
be contained in the static stress intensity factors for the 
bifurcated configuration given by H 33 (X, 0)K n . Indeed, 
Fig. 4 in ref. 12] shows that ratio 



H 33 {\,v'/c s 
H 33 (A, 0) 



(9) 



is very close to unity (up to ±5%) for all the values of 
A and v'/c s . Therefore, we conclude that even for this 
complex configuration the dynamic stress intensity factor 
admits an approximate decomposition in the form of a 
product of its static counterpart and a universal function 
of the local crack tip velocity, 



K' » y/l ~ v'/cs #33 (A, 0) K 



(10) 



This result suggests that a very good approximation for 
the dynamic stress intensity factor can be obtained by 
calculating the static stress intensity factor for the same 
instantaneous configuration and a knowledge of a univer- 
sal velocity function characteristic of the local symmetry 
conditions at the crack tip. 



The Decomposition Approximation for our 
Problem 



For advancing the problem posed in this paper we need 
to consider the stress field in the vicinity of three crack 
tips. In the vicinity of the tip of the macro crack we write 
(in the local polar coordinates (r M , M ) around that tip) 



£y(0 M) « M (t)) 



(11) 



where K M is the stress intensity factor that in principle 
depends on the positions and velocities of all the tips 
[i.e. L(t),l±(t),v u (t),v±(t)] and maybe other deriva- 
tives. Near the tips of the micro crack we write similarly 
(in local polar coordinates (r±,9±) around each tip) 



<rij(r±,0±,t) 



„ £y(0±,«±(t)) 
K ± 7x= 



(12) 



where K± are again the stress intensity factors 
that depend on all the time dependent functions 
L(t),£±(t),v M (t),v±(t) and maybe other derivatives. 

Our basic approximation is now motivated by the two 
examples Eqs. © and ifTUI) ; we assume that the dynamic 
stress intensity factors can be decomposed according to 



K 4 



fc I ( U+ )^(i(t),£ + (t),£_(t) 



(13) 



Here the universal function fc : (y) is the same function 
appearing in Eq. and all the stress intensity factors 
with superscript s refer to the solution of the static prob- 
lem with a frozen geometry which is given by the crack 
tip positions L(t),£±(t). On physical grounds we expect 
this approximation to be good when £/L — > 0, and to 
lose its validity as this ratio increases. The numerical 
simulations presented in Sect. Ivl lend a strong support 
to this expectation. 

The advantage of this approximation is that it leads 
to ordinary differential equations for the tip positions in 
much the same way that Eq. JBJ followed from Eq. ©, 



dL{t) 
dt 

d£-(t) 
dt 

d£+(t) 
dt 



c R 



cr 



cr 



ET 





e+(t),t-(t))]* 




ET 


[K'_(L(t), 


e+(t),t-(t))]* 




ET 



[K' + (L(t),t + (t),£-(t))f 



.(14) 



We turn now to the analysis of this set of equations and 
their consequences. 
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IV. SOLUTION OF THE MODEL 
A. The static problem 

A prerequisite to the solution of the set of equations 
(|14fl is the calculation of the static stress intensity fac- 
tors for a general configuration of two colinear cracks. 
We employ the available solution for two colinear cracks 
consisting of segments a < x < b and c < x < d with 
b < c under a remote mode I loading a°° . The a yy com- 
ponent of the stress tensor along the cracks line, outside 
the cracks, is given by 



a yy (x,0) 



2G(x) 



2x 2 - (a + b + c + d)x 



+ ab + cd — (d — b)(c — a) 



E(m) 



K(m) 



(15) 



Here 



G(x) 



y/ (x — a)(x — b)(x 

(d-c)(b-a) 

(d-b)(c-a) 



c)(x — d) 



(16) 



and E and K are the complete elliptic integrals of the 
first and second kind ^4[. The stress intensity factor at 
any one of the tips is obtained by taking the limit 



K i = lim y/2n(x - Xi)a yy (x,0), 



(17) 



where Xi is any one of positions of the tips. 

In order to adapt the general configuration to our 
macro crack and micro crack configuration we set a = 
— L, b = 0, c = d = £ + . Taking the limits in Eq. 
(|17f) . under the assumption L » <, we can extract the 
stress intensity factors at the three tips 



K, = <t°°a — 




(18) 



I- K(l -£-/£+) ^l +J 
E(l-i_/^ 



1 - 



K(l -£-/£+) J^l-£_/£ 



_ - 1 
1 



Note that the common pre-factor er°° ■J 'ttL /2 is just the 
stress intensity factor of the macro crack in the absence 
of the micro crack and serves here as the scale of the three 
stress intensity factors. In Fig. |3we present the three 
stress intensity factors as a function of A. In this example 
we kept the macro tip at L and the right micro tip at 
£+ fixed while £ was changed. We note that the stress 
intensity factor of the macro crack goes to the single crack 
result (unity in the reduced coordinates of Fig. |3J when 
£/A — » 0. Similarly, the stress intensity factor at £ + goes 
to unity when A — > 0, since also in that limit we remain 
with one crack. This last fact is not easily seen in Fig. 
13 since the upturn towards unity is very rapid, occurring 
just before coalescence. 
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FIG. 3: The normalized static stress intensity factors (SIFs) 
as a function of A. The normalization factor is <t°° U ttL/2, 
which is the stress intensity factor of the macro crack in the 
absence of the micro crack. We fixed £+ = 1 and varied 
Note that the stress intensity factors obey K+ < K- < K M 
and that K M — > <t°° */ tyL/2 as the ratio £/A decreases, as 
expected. 
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FIG. 4: The three crack tip velocities for an interaction event 
when a macro crack travelling at an initial velocity v M = 
0.62c_r interacts with a colinear micro crack of length £ = 
5 positioned at A = 5. The figure shows the normalized 
velocities v/cr as a function of the normalized time cnt/£. 



B. The dynamic problem 

Using the static stress intensity factors Eqs. (|18|l m 
Eqs. I|14[l we can solve numerically for the dynamics of 
the three tip positions. An example of the ensuing dy- 
namics is exhibited in Fig. 21 We note that what is seen 
in this picture is typical to all the conditions that we have 
considered: the macro crack is first accelerated, then the 
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FIG. 5: The normalized experimental velocity v exp /cR as a 
function of the normalized time cut /I. The experimental ve- 
locity was calculated according to Eq. 119H with A c <C I. 



left tip of the micro crack meets the fracture criterion 
Eq. and accelerates towards the macro crack; after 
some time lag, the right tip meets the fracture criterion 
and starts to move and attains at coalescence a lower 
velocity than the original macro crack. 

To connect to velocity fluctuations observed in experi- 
ment we reinterpret the data in Fig. 0]as it would be seen 
by an observer. We are physically motivated by the fact 
that as a result of a finite measurement resolution, below 
a critical separation A c the macro crack tip and the outer 
tip of the micro crack are indistinguishable and the mea- 
sured velocity is a result of some averaging. Therefore, 
let us define the experimental velocity v exp as 

v = { V ™ for A>A c 

P ~ I K + v+)/2 for A< A c 

Figure |3] shows the experimental velocity v exp during an 
interaction event. The rise in velocity after the steep de- 
cline is the second of Eqs. (|19|l . The last branch in which 
the velocity returns to the pre-collision value is out of the 
scope of the present model and had been added by hand 
for the sake of illustration. The point to stress is that 
the dynamic interaction event generates a typical large 
and rapid velocity "fluctuation" . It should be noted that 
we do not consider here the effect of the nucleation of 
the micro crack on the macro crack velocity. Physically 
we expect this effect to produce a sudden deceleration 
of the macro crack prior to the effect shown in Fig. |3l 
this is expected since the energy supply to the crack tip 
region should be partitioned between the nucleation pro- 
cess of the micro crack and the fracture process of the 
macro crack. This effect will be taken into account in 
future work where our current model will be coupled to 
a reasonable nucleation theory. 



V. SIMULATION AL SUPPORT 

Since our reduction to ordinary differential equations 
rests on the assumption of the product structure for 
the dynamic stress intensity factors, we must test the 
quality of the approximation by numerical simulations. 
We employ lattice simulations as described bellow. 

A. Lattice simulations 

Lattice models [H E E3 E H HJ provide a 
convenient, concrete, and physically sensible method of 
realizing crack dynamics. The material is represented by 
a lattice of mass points connected by Hookean springs. 
Fracture is achieved when a spring exceeds a certain crit- 
ical extension. In certain special cases, analytic solutions 
for static and steadily moving cracks can be obtained. 
In general, the model can be easily simulated. A ma- 
jor advantage of this class of models is that the process 
zone is quite small, on the order of a few lattice spac- 
ings. Thus, already on scales of 50 or so lattice spacings, 
continuum dynamics is very well realized. It has been 
recently demonstrated |22| that the universality assump- 
tion underlying linear elastic fracture mechanics, namely 
that the instantaneous crack velocity is only a function of 
the stress intensity factor at that moment, is extremely 
well satisfied by the lattice dynamics. 

For our present purposes, we use the machinery devel- 
oped in [23] ■ We work with a square lattice with nearest- 
and next-nearest- neighbor bonds, with the ratio of the 
spring constants chosen to give isotropic elasticity. As we 
are interested in cracks that propagate along the midline, 
we allow only bonds that cross the midline to break. We 
start with a lattice under fixed-displacement loading at 
the top and bottom, with bonds broken according to the 
desired initial configuration, Fig. ^ We relax this lattice 
using a multi-grid approach to accelerate convergence. 
At this point, we manually break the bond at the end of 
the macro crack, and monitor the subsequent sequence 
of bond breakings. 

B. Results of lattice simulations 

In order to test the quality the product structure ap- 
proximation we should first show that the analytic equa- 
tion (jSJ) for the dynamics of a single macro crack describes 
correctly the corresponding dynamics in the lattice sim- 
ulations. Multiplying Eq. © by L we obtain 

L^ = aL + P, (20) 

CR 

where a is predicted to equal one. The parameter [3 re- 
lates the material properties and boundary conditions of 
the lattice experiment to those in the continuum model. 
We simulated a single crack propagating without any ad- 
ditional damage, and measured v m /cr as a function of L. 
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FIG. 6: The simulation data of Lv m /cr as a function of L 
is shown in circles and the linear fit for this data is shown in 
solid line. It is clear that the functional form suggested by 
Eq. I20H indeed describes the single crack dynamics in the 
lattice simulation. The fit is consistent with the assumptions 
since the slope a is very close to one. 
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FIG. 7: The simulation data of v m /cn as a function of L is 
shown in circles and the analytic approximation predictions 
are shown in solid line. In this simulation we used L ps 245 
(the length of the macro crack in the interaction region) and 
I = 10. The analytic approximation predicted the simulated 
data up to an error of 6.5%, even though the ratio L/£ was 
not very large. 



Next we fitted the data to Eq. (|2(J|) . Figure shows that 
the functional form given by Eq. I|20[l describes well the 
single crack dynamics in the lattice simulation. More- 
over, we found that awl; thus our procedure appears 
internally consistent. It should be noted that, notwith- 
standing the excellent agreement evidenced in Figure El 
there are a number of uncontrolled approximations at 
play here. First, the functional form in Eq. (|2(Jfl is based 
on the assumption of a constant fracture energy. In fact, 
the fracture energy for our theoretical "lattice" material 
has been calculated, and it is not constant. Second, the 
simulation employs constant displacement boundary con- 
ditions in a finite strip (though, in order to mimic infinite 
medium, we specialized for times that do not allow wave 
interactions with the outer boundaries), and the theory 
assumes fixed stress at infinity. The excellence of the fit 
despite all this is somewhat unexpected, and bears fur- 
ther study. 

To directly test our product structure approximation 
Eqs. we used the value of j3 obtained by the linear 
fit and performed lattice simulation where a macro crack 
interacts with a colinear micro crack. We encounter a dif- 
ficulty in the simulations since the micro crack tips were 
trapped even when the stress intensity factor exceeds the 
material threshold . This known phenomenon of lat- 
tice trapping is an artifact of the lattice structure; to 
overcome it we fixed the micro crack tips also in the ana- 
lytic calculation. An example of the comparison between 
the simulation data and the analytic approximation is 
shown in Fig. We found that our analytic approxima- 
tion agrees with the simulated data, with deviations that 
are typically small. The largest errors are smaller than 
about 6 — 7% even for L/Ik, 25. For larger ratios, which 



is the expected physical regime, we expect better approx- 
imations. Note that the fact that the analytic approxi- 
mation overestimates the velocity of the macro crack is 
expected on physical grounds. Our product structure ap- 
proximation relies on the fact that for short distances the 
information on the positions of the cracks tips, carried by 
elastic waves, flows almost instantly even if the typical 
cracks propagation velocities are of the order of cr. In 
reality it takes finite time for the stresses to reorganize 
themselves according to the new cracks tips positions and 
generally the energy release rate is lower than in our ap- 
proximation. 

The main conclusion of this section is that the prod- 
uct structure approximation Eqs. (|13f) gives very good 
predictions for the model problem studied here. We pro- 
pose to interpret this in a broader way, and to test the 
applicability of this approximation in other contexts. 



VI. SUMMARY AND CONCLUSIONS 

This study has two aims; on the one hand, we are in- 
terested in the velocity fluctuations seen in dynamical 
crack propagation, and proposed here that the linear- 
elastodynamics interactions with micro cracks may very 
well be responsible for them. We did not address the 
reasons for the existence of the micro cracks - these may 
be there ab-initio or get born by the high stresses near 
the tip of the advancing cracks. On the other hand, we 
are after simplified methods of analysis of crack propaga- 
tion in non-trivial environments. The technical difficulty 
of solving the full dynamical equations calls for approxi- 
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mate methods that work. With the motivation presented 
in Sect. 1 1 i 1 1 we demonstrated how the assumption of 
product structure for the stress intensity factors reduces 
the dynamics to a set of ordinary differential equations 
that are easily solved. The gratifying agreement with 
the lattice simulations emboldens us to propose this as 
an approach that may find applications in other contexts 
of interest. Only future work will help to strengthen or 
delineate the usefulness of this approach. 



To further connect the interaction model to experimen- 
tal observations one should couple our theory to a phys- 
ically motivated model that will determined the condi- 
tions for micro cracks nucleation. Such a model will po- 
tentially predict the appearance of distributed damage 
in the process zone and the interaction with the macro 
crack may be related to the roughness of crack surfaces 
and the quasi-periodicity of the velocity fluctuations. 
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